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The configurational and melting properties of large two- 
dimensional (2D) clusters of charged classical particles inter- 
acting with each other via the Coulomb potential are inves- 
tigated through the Monte Carlo simulation technique. The 
particles are confined by a harmonic potential. For a large 
number of particles in the cluster (N>150) the configuration 
is determined by two competing effects, namely in the center 
a hexagonal lattice is formed, which is the groundstate for 
an infinite 2D system, and the confinement which imposes its 
circular symmetry on the outer edge. As a result a hexago- 
nal Wigner lattice is formed in the central area while at the 
border of the cluster the particles are arranged in rings. In 
the transition region defects appear as dislocations and discli- 
nations at the six corners of the hexagonal-shaped inner do- 
main. Many different arrangements and type of defects are 
possible as metastable configurations with a slightly higher 
energy. The particles motion is found to be strongly related 
to the topological structure. Our results clearly show that 
the melting of the clusters starts near the geometry induced 
defects, and that three different melting temperatures can be 
defined corresponding to the melting of different regions in 
the cluster. 

PACS numbers: 45.05.+X, 61.46.+W, 73.22.-f 



I. INTRODUCTION 

Recently, there has been considerable theoretical and 
experimental progress in the study of mesoscopic sys- 
tems consisting of a finite number of charged particles 
which are confined into an artificial circular symmetric 
potential. In 1934, Wigner suggested that a liquid to 
solid phase transition should occur in a three-dimensional 
(3D) Fermi system at low densities JjJ. Typical experi- 
mental model systems for the study of this system are 
electrons on the surface of liquid helium , electrons in 
quantum dots H, colloidal suspensions [H and in con- 
fined plasma crystals [^|. On the other hand, various 
similar systems, like the vortex clusters in an isotropic 
superfluid 0, vortices in superfluid He 4 [J?], vortices in 
a Bose-Einstein condensate stirred with a laser beam || 
and in mesoscopic superconducting disks |^| have many 
common features with those of 2D charged particles. Col- 
loidal particles dissolved in water |lfj|] are another exam- 
ple of an experimental system where classical particles 
exhibit Wigner crystallization. Recently, macroscopic 
2D Wigner islands, consisting of charged metallic balls 
above a plane conductor were studied and ground state, 
metastable states and saddle point configurations were 
found experimentally JXTt] - 



Such a system with a finite number of particles, ini- 
tially studied by Thomson as a classical model for the 
atom |l^,|l3|, has been extensively studied during the 
past few years. For a small number of particles (typi- 
cally N < 100) they are arranged in rings |l4|~[l7j and 
a Mendeleev-type of table was constructed in Rcf. pl| 
which gives the distribution of those particles over the dif- 
ferent rings. Moreover, the configurations of the ground 
state, the metastable states and saddle point states were 
obtained, from which the transition path and the geomet- 
ric properties of the energy landscape were given in Ref. 
fl8f . The spectral properties of the ground state config- 
urations were presented in Ref. |16| and generalized to 
screened Coulomb 19 2^] and logarithmic Jf|^0| inter- 
particle interactions. The excitation of normal modes of 
2D Coulomb clusters in laboratory complex plasmas were 
recently observed [ Ml . 

The melting properties of this system have been stud- 
ied experimentally (l(]j2§] and by Monte Carlo(MC) 
p3| , p4j and molecular dynamics [^5|,^6| simulations. In a 
hard wall confined system with short-range inter-particle 
interaction, the melting behavior was found even more 
interesting. Reentrant melting of 2D colloidal clusters in 
a hard wall potential was obtained in both experimental 
|Tof and theoretical work [^5| . 

The defect structure in crystals is of paramount im- 
portance for the stability and the strength of mate- 
rials. Topological defects in Wigner crystals [^7 28 
and its effect on particle melting were investigated in 
Refs. [ p2| , p3p^ , |29| ] . Thermal defect mediated melting was 
proposed as the microscopic mechanism for melting in an 
infinite 2D triangular Wigner crystal. It is well known 
that the KTHNY scenario(Kosterlitz-Thouless-Halperin- 
Nelson- Young) describes 2D melting as a defect-mediated 
phenomenon and melting occurs in two stages through 
two continuous phase transitions |30|. A clear case of 
two-stage melting was observed in a film of paramagnetic 
colloidal particles ^2| . 

In this paper we study topological defects which are 
induced by the confinement potential, i.e. which are a 
result of the finite size of the system. Next we investigate 
how these defects influence the melting of the mesoscopic 
2D island. The present paper is organized as follows. In 
Sec. II, we describe the model system and the numerical 
approach. Sec. Ill is devoted to the structural proper- 
ties of the topological defects at zero temperature. In 
Sec. IV, we discuss the eigenmode spectrum for these 
large clusters. The discussion on the non-homogeneous 
melting is presented in Sec. V. Our conclusions are given 
in Sec. VI. 
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II. NUMERICAL APPROACH 

The model system was defined in Ref. jl5| and the 
Hamiltonian for such a system is given by 
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m*oj 2 r 2 is taken 



The confinement potential V(r) = i 
circular symmetric and parabolic, where m* is the effec- 
tive mass of the particles, q is the particle charge, luq 
is the radial confinement frequency and e is the dielec- 
tric constant of the medium the particles are moving in. 
To exhibit the scaling of the system, we introduce the 
characteristic scales in the problem: tq — (2q 2 /meuj 2 .) 1 ^ 
for the length, Eq = (muj 2 q 4 /2e 2 ) 1 ^ 3 for the energy and 
To = (mujQq 4 /2e 2 ) 1 / 3 kg 1 for temperature. After the scal- 
ing transformations (r — » r/r ,E — > E/E ,T — > T/Tq), 
the Hamiltonian can be rewritten in a simple dimension- 
less form as 
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which only depends on the number of particles N. The 
numerical values for the parameters Wo, ?"o> Eq, Tq for 
some typical experimental systems were given in Ref. [ p"5[ . 

The MC simulation technique |3l|] is relatively simple 
and rapidly convergent and it provides a reliable esti- 
mation of the total energy of the system in cases when 
relatively small number of Metropolis steps is sufficient. 
However, the accuracy of this method in calculating the 
explicit states is poor for systems with a large number 
of particles, which have significantly more metastable 
states. To circumvent this problem we employed the 
Newton optimization technique which was outlined and 
compared with the standard Monte Carlo technique in 
Ref. p6| . The structure and potential energy of the 
system at T ^ are found by the standard Metropo- 
lis algorithm |31| in which at some temperature the next 
simulation state of the system is obtained by a random 
displacement of one of the particles. We allow the system 
to approach its equilibrium state at some temperature T, 
after executing 10 4 -5 x 10 5 "MC steps" . Each MC step is 
formed by a random displacement of all particles. If the 
new configuration has a smaller energy it is accepted, but 
if the new energy is larger the configuration is accepted 
with probability 5 < exp(— AE/T), where 8 is a random 
number between and 1 and AE is the increment in the 
energy [pi. 



III. TOPOLOGICAL DEFECTS 

It is well known that the hexagonal lattice is the most 
energetically favored structure for classical point charges 



in a two-dimensional infinite plane at low temperature 
fS2f . For a system consisting of a finite number of re- 
pelling particles restricted to 2D, which are held together 
by a circular harmonic potential, the cluster patterns 
are determined by the need to balance the tendency to 
form a triangular lattice against the formation of a com- 
pact circular shape. The configuration is determined by 
these two competing effects, namely circular symmetry 
and triangular structure (Wigner lattice). This competi- 
tion leads to intrinsic defects in the 2D circular Coulomb 
cluster which are geometry (of the confinement poten- 
tial) induced defects. This ground state is not a defect 
free system. The symmetry breaking is due to the pack- 
ing of the triangular lattice into a region with a circular 
boundary. A hexagonal lattice which is cut by a cir- 
cle without the introduction of any defect has an energy 
E = 56.0499i?o which is larger than the ground state 
energy E = 55. 9044E for N = 291 particles. 

In the first part of this paper, we investigate the form 
and position of the defects in large clusters. Therefore we 
make use of the Voronoi construction |33| . The Voronoi 
construction of a collection of particles consists of a parti- 
tion of space into cells. Each cell consists of those points 
which are closer to the given particular particle than to 
any other particles. Examples of Voronoi constructions 
are shown in Fig. ^ where the ground state configuration 
for N = 291,300,400 and 500 are shown. One can see 
that there are two kinds of defects, i.e. dislocations and 
disclinations. Disclinations are orientational defects with 
five (indicated by '-') or seven (indicated by '+') fold co- 
ordination number (the number of sides of the polygon 
around the particles is nothing else then the coordination 
number). A dislocation is a pair of two disclinations con- 
sisting of a defect with 5-fold (-) and a defect with 7-fold 
(+) coordination number. In the latter case the order- 
ing at long distances is not disrupted and consequently 
such a bound pair has a much lower energy ifjof . The 
total number of 5-fold A_ and 7-fold N + disclinations 
depends on the particular configuration. The number of 
disclinations in this system is determined by Euler's the- 
orem and can't be changed, so the net topological charge 
N- — N + is always equal to six as was already demon- 
strated in Refs. The reason is that every '-' defect 
can bend the lattice clockwise over 7r/3 from a straight 
lattice and thus six '-' defects can bend a straight line 
into a circle. Dislocations will appear when it decreases 
the energy of the system. From Fig. 1 it is apparent that 
this is more so for larger clusters. 

In Refs. 20 2(|, the defects in clusters with a logarith- 
mic inter-particle interaction were studied. We want to 
stress that their way of visualizing the defects is differ- 
ent: nearest neighbours are connected by a line, without 
making crossings. However, this does not lead to a unique 
picture: the total number of 5-fold, JV_, and 7- fold, N+, 
disclinations can vary in the same configuration, only the 
net topological charge A_ — N + is always equal to six. 

In these large clusters, the defects are located on a 
hexagon, i.e. they form a hexagonal structure. As can 
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be seen in Fig. [fj the defects are approximately situated 
at the six corners of a hexagon, each corner with a net 
topological charge of Notice from Fig. that a single 
5-fold disclination can appear, but never a single 7-fold 
disclination. For the large strain energy around '— 1' topo- 
logical charge, some dipole defects(i.e. dislocation with 
'— 1' and '+1' defects) will be generated to shield the '- 
1' topological charge. Those shielding dipole defects do 
not change the topology of the system. A clearer ex- 
ample is shown in Fig. || for the N = 291 groundstate 
configuration. We considered the N — 291 system as it 
minimizes the number of defects. The reason is that for 
this particle number the configuration has 42 particles 
in the outer ring, which is a multiple of the topological 
charge. There are three rings at the border with an equal 
number (N = 42) of particles (the ID Wigner lattice), 
the central hexagonal structure (the 2D Wigner lattice) 
and the defects indicated by triangles (A) and squares 
(□) are situated: (i) around the six corners of a hexagon 
and (ii) in the transition region between the outer rings 
and the central hexagon. 

It should be noted that the search for the global min- 
imum configuration is a difficult problem for large sys- 
tems because of the existence of a large number of local- 
minimum configurations, with energy very close to the 
global minimum. Thus one is never 100% sure to have 
found the real ground state. Therefore, we investigated 
the different metastable states. In an experiment, those 
metastable states may be reached by thermal excitation 
if the energy barrier between them and the ground state 
is comparable to or smaller than ksT. The saddle points 
between those metastable states were investigated in an 
earlier paper @ for N < 20. In Fig. | the energy and 
the total number of defects of different metastable states 
are shown for N = 300. The results for the different 
metastable configurations are ordered with increasing en- 
ergy. Note from Fig. 0(b) that on average the total num- 
ber of defects increases with energy, but it shows strong 
local variations. Only an even number of defects are ob- 
tained, because the net topological charge is always six, 
and the dipole defects (i.e. one dislocation with '-1' and 
'+1' defects) always appear in pairs. Also the hexago- 
nal position of the defects disappears (see Fig. 0(1)) and 
more free dislocations are found. These defects move 
from the transition region to the border (see Fig. 0(2)) 
or to the central region (see Fig. 0(3)). For confi gura- 
tions with higher energy, the defects arrange themselves 
in long chains, i.e. dislocation lines. On average the con- 
figurations with defects at the border have a lower energy 
than those with defects in the center. 

We also investigated whether or not it is possible to 
have a configuration with only six 5-fold disclinations 
and no other defects (like for example for the N = 85 
configuration with 24 particles at the outermost ring 
p4j). Therefore, we started our MC procedure with a 
perfect hexagonal structure without any defect and then 
allowed it to relax to its energy minimum. We did this 
for N = 281 up to 295 particles, because we noted that 



for these particle numbers the configuration has about 42 
particles in the outer ring, which is a multiple of six, i.e. 
the net topological charge. Only in such a case one can 
have the situation in which just six 5-fold disclinations 
are present. We found that our result (from N = 281 up 
to 295 particles) never converges to a configuration with 
only six 5-fold disclinations. However, this procedure in- 
deed favorably relaxes to configurations with 42 particles 
in the outer ring, often resulting in a configuration which 
has less total number of defects than the corresponding 
groundstate. 

IV. THE EIGENMODE SPECTRUM 

The effect of the geometry induced defects on the 
eigenfrequencies (i.e. the eigenmode spectrum) were also 
investigated for these large clusters. In this system, it is 
well known that there are three eigenfrequencies which 
are independent of N |l6): uj = 0, \/2 and \f&, which cor- 
respond to the rotation of the system as a whole, the cen- 
ter of mass mode and the breathing mode, respectively. 
The above modes were recently observed experimentally 
|^lf . The smallest frequency no longer correspond to in- 
tershell rotation as in small clusters but to the excita- 
tion of a vortex/antivortex pair, of which a typical mode 
is shown in Fig. 0(a). Slightly larger excitation energies 
may consist of multiples of such pairs (see Fig. 0(b)). 
Modes with higher eigenfrequencies often show a hexag- 
onal structure similar to the ordering of the defects. The 
motion can be concentrated around (see Fig. [|(c)) or be- 
tween the defects (see Fig. 0(d)). The local modes can be 
found at the six corners of the hexagon where the defects 
are exactly situated (see Fig. 0(e)). The modes in which 
the inner particles have larger amplitudes than the outer 
particles have the largest eigenfrequencies (see Fig. 0(f)). 

The lowest eigenfrequencies of the excitation spectrum 
corresponding to the ground-state configuration of the 
system is shown in Fig. 0, as function of the number of 
particles for N ranging from 281 to 307. The labels in 
Fig. denote the total number of defects present in the 
ground state. Notice that only an even number of defects 
are obtained as explained before. On average, configura- 
tions with a large number of defects have a smaller lowest 
eigenfrequency and are thus less stable, and vice versa. 

In this 2D lattice, all behaviors of the cluster modes can 
be classified as shear-like or compression-like modes. In 
order to characterize the compressional and shear parts 
of these eigenmodes, we calculated respectively the di- 
vergence V- v and the vorticity (Vx v) z of the velocity 
field. To calculate the velocity field, we interpolated the 
displacements of Fig. 4 on a 100 x 100 grid (thus ne- 
glecting the constant eigenfrequency). The "divergence 
and vorticity maps" were then calculated at every point 
of this matrix. Notice that pure shear or compressional 
modes do not exist in the circular boundary of finite clus- 
ter. Fig. 6(b) shows the vorticity and thus displays the 
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shear part of the eigenmode of Fig. 4(b). The two vor- 
tex/antivortex pairs are clearly seen. The "divergence 
map" for this eigenmode is practically zero everywhere, 
as there is no compressional part (Fig. 6(a)). This is not 
the case for the eigenmodes (c) and (d) of Fig. 4. Figs. 
6(c) and (e) show the "divergence maps" for both eigen- 
modes, in which the compression and rarefaction can be 
clearly seen. Both cases show no shear part (Figs. 6(d) 
and (f)). Figs. 6(c) and 6(e) (see also the 3D plots at the 
bottom of Fig. 6) exhibit clearly dipole type of compres- 
sional oscillations between (Fig. 6(c)) and at (Fig. 6(e)) 
the defect regions. 

V. NON-HOMOGENEOUS MELTING 

Understanding the microscopic mechanism of melting 
has intrigued scientists since the late nineteenth century. 
Special interest has been devoted to 2D melting [p5[ . 
Most works address infinite systems consisting of a single 
layer. However, whether melting of a 2D crystal is a first 
order transition and proceeds discontinuously or is a con- 
tinuous transition in which the crystal first transits into 
a hexatic phase retaining quasi-long-range orientational 
order and then melts into an isotropic fluid, is still an 
open question and a controversial issue. 

In the present work we consider a finite 2D system 
where we take N — 291 for our numerical simulation. 
Here we present a calculation of the meltin g ph ase dia- 
gram by performing MC simulations. In Ref. |2q1 molecu- 
lar dynamics was used to investigate the melting of a clus- 
ter of particles interacting through a logarithmic interac- 
tion. As compared to our Coulomb interaction where the 
geometry induced defects are situated in the 3*^ and 4 th 
outer shells (i.e. the transition region) and around the 
six corners of the "defect" hexagon, in the logarithmic 
interacting system p(| those defect are mainly situated 
in the outer two shells. In Ref. [2q] , the number and 
type of defects were studied as function of the noise (i.e. 
temperature). Here we will use several different crite- 
ria such as the total energy, the radial dependent mean 
square displacement, the bond- angular order factor and 
the angular square deviation to characterize the melting 
behavior of the cluster. 

There are several different criteria that can be used to 
find the melting temperature. In order to determine the 
melting transition point, we calculated the potential en- 
ergy of the system as a function of temperature (see Fig. 
^) . In the crystalline state the potential energy of the sys- 
tem increases almost linearly with temperature and then 
after the critical temperature is reached (T/Tq = 0.01 
for N — 291), it increases more steeply as shown in Fig. 
g. This is a signature of melting and is related to the 
unbinding of dislocation pairs. The dotted assurgent line 
in Fig. 7 indicates the linear temperature dependence of 
the potential energy for low temperatures before melting. 
In the upper inset Fig. fz](a), we plot AE which is the dif- 



ference between the numerical obtained energy and the 
linear T-behavior. After the melting point, AE increases 
sup- linear. 

The lower inset, Fig. ||(b), shows the averaged num- 
ber of defects as function of temperature T/Tq. The 
number of defects were obtained as follows. We con- 
sidered 40 configurations for every temperature T/Tq. 
Every 500 MC steps a new configuration was obtained. 
For all these configurations the number of defects were 
counted. Finally we averaged over the 40 configurations, 
which is the reason why the number of defects can be 
non-integer. With increasing temperature, the system 
generates more and more defects and after the melting 
point the defect number grows very fast. Notice that two 
clear critical temperatures emerge from this figure at the 
crossing points of the dotted lines, i.e. T/Tq = 0.01 and 
T/T a = 0.014. 

In Fig. 8 we plot typical particle trajectories for dif- 
ferent temperatures which shows that the melting of this 
system is very complex and non-homogeneous. It clearly 
indicates that the melting starts around the six corners of 
the hexagon which are exactly the defect regions. With 
increasing temperature, the particles in the defect region 
start to move radially and destroy order locally. With 
further increase of temperature the total system com- 
pletely melts and the order is destroyed. 

In order to better describe the spatial dependence of 
the melting process in our system, we separate the con- 
figuration into three regions as shown in Fig. ^. Region 

I (dark grey colored hexagonal area) is comprised of the 
defect-free hexagonal center; region II is a transition re- 
gion with the defects (light grey colored area), and region 
III consists of the outermost two rings. For the case of 
N = 291 particles region I consists of 91 particles, region 

II consists of 116 particles and region III of 84 particles. 
We calculate for each region the mean square displace- 
ment (ufj), which was introduced in Ref. [fl5|| , 

1 N 

("«) = ^E(^-^)) 2 )/« 2 ' ( 3 ) 

with a = 2R/y r N the average distance between the par- 
ticles. Fig. shows the (wfj) as a function of the re- 
duced temperature T/Tq for the three different regions. 
At low temperatures the particles exhibit harmonic oscil- 
lations around their T = equilibrium position, and the 
oscillation amplitude increases linearly and slowly with 
temperature: the particles are well localized and display 
still an ordered structure. This linear dependence is ac- 
centuated by the thin straight lines in Fig. 9. Here, we 
already notice that the amplitude of the local particle 
thermal vibrations in these different regions are differ- 
ent. The amplitude is largest at the defect region and 
lowest in the center of the cluster. Melting occurs when 
(u^) increases very sharply with T. Because of the fi- 
nite number of particles one has rather a melting region, 
instead of a well-defined melting temperature. After the 
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melting ''point", the particles exhibit liquid-like behav- 
ior. Fig. |9j exhibits three different melting temperatures 
corresponding to the three different regions. Firstly re- 
gion II, i.e. the transition region containing the defects, 
starts to melt, then the outermost two rings melt, and fi- 
nally the hexagonal region melts. Following Ref. P^| , we 
can "define" a melting temperature at the point where 



I) 



0.10 which results into the melting temperatures 



T meU /T ~ 0.0115,0.0125 and 0.0138 for the defect re- 
gion, the outer rings and the center region, respectively. 

In order to investigate the melting in the defect region 
in further detail we consider two new small regions as 
shown in the inset of Fig. [To]. One region is around the 
defect, the other doesn't contain a defect and is situated 
between two defect regions. For N — 291 each of the 
two regions contains respectively eight and seven parti- 
cles. In Fig. |h], the («#) of these two different regions 
show a different melting temperature: the melting clearly 
starts first around the defect as expected. The particle 
motion is strongly influenced by the topological defects, 
i.e. the particles in the defect regions are less well in- 
terlocked and have a larger diffusion constant than the 
undistorted lattice regions and their thermal motions are 
easier to be excited [^9|. Notice, that for the two sep- 
arate regions a much sharper melting behavior is found 
than for the intermediate region as a whole (see Fig. ||). 
The reason of course is that in Fig. ^ one averages over 
defect and defect-free regions. The criterion (u\\ » 0.10 \ u ai) 
results into T me i t /T ~ 0.0118,0.0138 for the defect and 
the defect-free regions, respectively. These two melting 
temperatures are very similar to the melting temperature 
of the transition region and the hexagonal region of Fig. 
9. 

The third independent parameter is the bond- 
orientational correlation function. This quantity deter- 
mines the type of melting transition and the melting 
point for an infinite system. Our finite system is too 
small in order to have a reliable analysis of the asymp- 
totic decay of the density correlation function. There- 
fore, we calculate the bond-angular order factor which 
was originally presented in Ref. [37[| , but following Ref. 
|E3i we modified it into, 



1 N 1 



N^N nb 



(4) 



This quantity is calculated only for region I which ex- 
hibits a hexagon structure, where j means the N n b near- 
est neighbors of particle i, for ideal hexagonal lattice 
N n b = 6, where 6j n is the angle between some fixed 
axis and the vector which connects the j th particle and 
its nearest n th neighbor. 

For a perfect hexagonal system Gq = 1. In our system 
for N — 291, the initial value of Gq is 0.96, which means 
that the structure in region I is almost perfect hexagonal. 
Our numerical results (see open dots in Fig. ||) show 
that Ge decreases linearly with increasing temperature. 
When Gq is around 0.6, it more rapidly drops to zero with 



increasing temperature. G§ should be zero for the liquid 
state. This can be compared with the infinite system 
where a universal melting criterion was found in Ref. p3| : 
melting occurs when the bond-angle correlation factor 
becomes Ge ~ 0.45, which was found to be independent 
of the functional form of the interparticle interaction. For 
our system the value Ge ~ 0.45 is probably not correct 
because in our finite system Gq does not drop to zero at 
T m eit, but is smeared out around T me it- Therefore, the 
midpoint G§ w 0.45/2 w 0.225 is expected to describe 
better the melting temperature. This leads to T me i t /T ~ 
0.0136 which is similar to the result T me /t/To ~ 0.0138 
obtained from the radial displacement criterion. 

In contrast to bulk systems, the melting scenario of 
small laterally confined 2D systems was found earlier (lj| 
to be a two step process. Upon increasing the tempera- 
ture, first intershell rotation becomes possible where ori- 
entational order between adjacent shells is lost while re- 
taining their internal order and the shell structure. At 
even higher temperatures, the growth of thermal fluctu- 
ations leads to radial diffusion between the shells, which 
finally destroys the positional order. To characterize the 
relative angular intrashell and the relative angular inter- 
shell, we use the functions as defined in Ref. fL5| . The 
relative angular intrashell square deviation 



1 



E 

i=l 



[(to 



Vii) ) - {<Pi - <Pn) /Vo, (5) 



and the relative angular intershell square deviation 



where i\ indicates the nearest particle from the same 
shell, while 12 refers to the nearest-neighbor shell, (fo = 
2ir/Nft, where the number in the outermost two rings 
Nr is the same and equals 42 for our N = 291 system. 
Only the two outermost rings have a clear shell structure. 
Both two outer rings are strongly interlocked which is a 
consequence of the ID Wigner lattice arrangement of the 
two rings. From the inset of Fig. [n], one can see that 
the inner ring will melt before the outermost ring. We 
find that the result for (u^i) of the inner ring is almost 
the same as (ti^ 2 ) which is the relative angular intershell 
square deviation. It means that when the inner ring loses 
its order, the relative order is lost simultaneously. The 
outermost ring can still keep its order and it will melt at 
even higher temperature. Comparing this with Fig. 9, 
we see that the radial and angular displacements start to 
increase rapidly at approximately the same temperature. 
Thus for large clusters intershell rotation will not occur 
below the melting temperature, but appears at the same 
temperature when the radial displacements increase. 
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VI. CONCLUSION 

The configurational and melting properties of large 
two-dimensional clusters of charged classical particles in- 
teracting with each other via the Coulomb potential were 
investigated through the Monte Carlo simulation tech- 
nique. For the ground state configuration, a hexagonal 
Wigner lattice is formed in the central area while at the 
border of the cluster the particles are arranged in rings. 
In the transition region between them defects appear as 
groups of dislocations and disclinations at the six cor- 
ners of the hexagonal-shaped inner domain. Many dif- 
ferent arrangements and types of defects are possible as 
metastable configurations with a slightly higher energy. 
The particles motion is found to be strongly related to the 
local topological structure. Our results clearly show that 
the melting of the clusters starts near the geometry in- 
duced defects, and that three melting temperatures can 
be obtained: T melt /T ~ 0.0115,0.0125 and 0.0138 for 
the defect region, the outer rings and the center region, 
respectively. These values are for the N = 291 cluster. 
Taking a different value for N does not lead to any qual- 
itative differences, it only influences slightly the values 
for the three melting temperatures. 
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FIG. 1. The ground state configurations for N=291, 300, 
400, 500 particles. The Voronoi structure is shown and the 
defects (i.e. disclinations) are indicated by '+' for a 7-fold 
and by '-' for a 5-fold coordination number. 

FIG. 2. The ground state configuration for TV = 291. The 
dots give the postion of the particles. Three regions are found: 
I (dark grey colored hexagonal area) is comprised of the de- 
fect-free hexagonal lattice; II is a transition region with the 
defects (light grey colored area), and III consists of the outer- 
most two rings. Obviously, there are three rings at the border, 
inner region forms the hexagonal lattice, and the defects are 
situated in the transition region, and they are at the exact six 
corners of the hexagon. The '+1' and '-1' topological defects 
are represented by the open squares and triangles, repectively. 

FIG. 3. The energy E/Eo (a) and the total number of de- 
fects (b) of different metastable states are shown for N = 300. 
Three typical defect configurations with different energy are 
shown at the right side of the figure. 

FIG. 4. Vector plot of the eigenvectors for the cluster with 
N = 291 particles for six different values of the mode number 
K. 

FIG. 5. Excitation spectrum of normal modes as a function 
of the number of particles in the cluster. The frequency is in 
units of u>'= wo/2 1 ' 2 . The numbers in the figure indicate the 
number of defects found in the ground state of the different 
clusters. 



FIG. 6. Gray-scale contour maps of the vorticity (Vx v) z 
and the divergence V- v of the velocity field of three different 
eigenmodes. A corresponding 3D plot is shown for those maps 
which exhibit a clear structure. 

FIG. 7. The potential energy ( E / E Q ) of the 2D Coulomb 
cluster as a function of temperature T /To for N = 291. The 
insets show AE = E — En ne (a), and the number of defects 
(b) as a function of temperature T/Tq. 



FIG. 8. Typical snapshots of particle trajectories for differ- 
ent temperatures T/T Q for N = 291. 

FIG. 9. The mean square displacements as function of the 
temperature T/To for the three regions defined in Fig. 2. 
The open symbols are the results for the correlation function 
refered to the right scale, for the inner hexagonal region. 
The linear dependence at low temperature is accentuated by 
the thin straight lines. The dotted curves are guides to the 
eye. 



FIG. 10. The mean square displacements as function of 
the temperature T/To for the small defect-free (open sym- 
bols) and defect regions (solid symbols) in the intermediate 
region as indicated by the circular areas in the inset. The 
thin straight lines show the low temperature linear depen- 
dence. The dotted curves are guides to the eye. 

FIG. 11. The relative angular intrashell square deviation 
(liai) an d relative intershell square deviation (u 2 a2 ) of the 
outermost two rings as a function of temperature for N = 291. 
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